library(tidyverse)
## -- Attaching packages ------------------------------ tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.1
## v tidyr   1.1.1     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## -- Conflicts --------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(getlandsat)
library(mapview)
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows

Question 1:

bb=read_csv("../data/uscities.csv") %>%
  filter(city=="Palo") %>%
  st_as_sf(coords= c("lng","lat"),crs=4326) %>%
  st_transform(5070) %>%
  st_buffer(5000) %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_as_sf()
## Parsed with column specification:
## cols(
##   city = col_character(),
##   city_ascii = col_character(),
##   state_id = col_character(),
##   state_name = col_character(),
##   county_fips = col_double(),
##   county_name = col_character(),
##   county_fips_all = col_character(),
##   county_name_all = col_character(),
##   lat = col_double(),
##   lng = col_double(),
##   population = col_double(),
##   density = col_double(),
##   source = col_character(),
##   military = col_logical(),
##   incorporated = col_logical(),
##   timezone = col_character(),
##   ranking = col_double(),
##   zips = col_character(),
##   id = col_double()
## )
#mapview(bb)

Question 2:

bbwgs= bb %>% st_transform(4326)

bb=st_bbox(bbwgs)


scene=lsat_scenes() %>%
  filter(min_lat<= bb$ymin, max_lat >= bb$ymax) %>%
  filter(min_lon<= bb$xmin, max_lon >= bb$xmax) %>%
  filter(as.Date(acquisitionDate) == as.Date("2016-09-26"))

write.csv(scene,file="../data/palo-flood.csv")


###

meta= read_csv("../data/palo-flood.csv")

files= lsat_scene_files(meta$download_url) %>%
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse= "|"), file)) %>%
  arrange(file) %>%
  pull(file)

st= sapply(files,lsat_image)

s=stack(st) %>% setNames(c(paste0("band", 1:6)))

#The dimensions of my stacked image are: 7811 (nrow)/ 7681 (ncol)/ 59996291 (ncell)/ 6 (nlayers)

#The CRS is: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 

#The cell resolution is: 30, 30  (x, y)

cropper= bbwgs %>%  st_transform(crs(s))

r=crop(s, cropper)

#The dimensions of my cropped image stacked are: 340 (nrow)/ 346 (ncol)/ 117640 (ncell)/ 6 (nlayers)

#The CRS is: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 

#The cell resolution is: 30, 30  (x, y)

Question 3:

par(mfrow=c(1,2))
#True Color
plotRGB(r, r=4 ,g=3 ,b=2, stretch= "lin")
plotRGB(r, r=4 ,g=3 ,b=2, stretch= "hist")

#Color Infared
plotRGB(r, r=5 ,g=4 ,b=3, stretch= "lin")
plotRGB(r, r=5 ,g=4 ,b=3, stretch= "hist")

#False Color Water Focus
plotRGB(r, r=5 ,g=6 ,b=4, stretch= "lin")
plotRGB(r, r=5 ,g=6 ,b=4, stretch= "hist")

#My Choice
plotRGB(r, r=1 ,g=3 ,b=6, stretch= "lin")
plotRGB(r, r=1 ,g=3 ,b=6, stretch= "hist")

#Color stretching determines how the values are represented on the image.

Question 4:

#Step 1

palette= colorRampPalette(c("blue","white","red"))
ndvi=(r$band5 - r$band4) / (r$band5 + r$band4)
ndwi=(r$band3 - r$band5) / (r$band3 + r$band5)
mndwi=(r$band3 - r$band6) / (r$band3 + r$band6)
wri=(r$band3 + r$band4) / (r$band5 + r$band6)
swi= 1 / sqrt(r$band2 - r$band6)
## Warning in sqrt(getValues(x)): NaNs produced
plot(ndvi,col=palette(256), legend=FALSE)

plot(ndwi,col=palette(256), legend=FALSE)

plot(mndwi,col=palette(256), legend=FALSE)

plot(wri,col=palette(256), legend=FALSE)

plot(swi,col=palette(256), legend=FALSE)

#Describe the 5 images. How are they simular and where do they deviate?

#All five images give us the results of different index methods which we first had to state. The normalized methods in short minimizes the effects of illumination and enhance spectral features that aren't originally visible. Our water threshold which is given show us areas that both retain and don't retain water & for NDVI index, numbers less than 0 = water shown in blue and numbers greater than 0 = vegetation shown in red. For both NDWI and MNDWI Index's, numbers greater than 0 = water shown in red and numbers less than 0 = vegetation without water presence. For WRI Index, it shows the water ratio by displaying water presence in cells with numbers greater than 1 in red. The same can be said for the SWI Index since it displays water presence in cells with values greater than 5.  


#Step 2

thresholding1= function(x){ifelse(x <= 0,1,NA)}
thresholding2= function(x){ifelse(x >= 0,1,NA)} 
thresholding3= function(x){ifelse(x >= 0,1,NA)} 
thresholding4= function(x){ifelse(x >= 1,1,NA)}
thresholding5= function(x){ifelse(x <= 5,1,NA)}

flood1= calc(ndvi,thresholding1)
flood2= calc(ndwi,thresholding2) 
flood3= calc(mndwi,thresholding3)
flood4= calc(wri,thresholding4)
flood5= calc(swi,thresholding5)

plot(flood1, col= "blue", legend=FALSE)

plot(flood2, col= "blue", legend=FALSE)

plot(flood3, col= "blue", legend=FALSE)

plot(flood4, col= "blue", legend=FALSE)

plot(flood5, col= "blue", legend=FALSE)

#mapview(flood1)
changing NA’s to 0s for Q5
thresholdings1= function(x){ifelse(x <= 0,1,0)}
thresholdings2= function(x){ifelse(x >= 0,1,0)} 
thresholdings3= function(x){ifelse(x >= 0,1,0)} 
thresholdings4= function(x){ifelse(x >= 1,1,0)}
thresholdings5= function(x){ifelse(x <= 5,1,0)}
floods1= calc(ndvi,thresholdings1)
floods2= calc(ndwi,thresholdings2) 
floods3= calc(mndwi,thresholdings3)
floods4= calc(wri,thresholdings4)
floods5= calc(swi,thresholdings5)

Question 5:

  #1. Split raster data/structure into just data
  #2. Calculate KMeans clustering on data
  #3. Add clustered data to original raster structure

set.seed(1)

#1. Splited Raster and get data to cluster
r_extract = getValues(r)
r_nna = na.omit(r_extract) #No NAs, no need to omit
#data was extracted band by band (there are 6 columns)

#2. Calcualting the KMeans Cluster
kmeans_r = kmeans(r_extract, 12)
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 5882000)
fviz_cluster(kmeans_r, geom="point", data = r_extract)

#3. Applying the cluster to raster structure
kmeans_raster = floods2
values(kmeans_raster) = kmeans_r$cluster
plot(kmeans_raster, col = viridis::plasma(12))

#Step 5.3

com_table = table(values(floods2), values(kmeans_raster))

which.max(com_table[2,])
## 12 
## 12
rast_thresh = function(x){ifelse(x != which.max(com_table[2,]),0,1)}
kmeans_floodmask = calc(kmeans_raster, rast_thresh)

floodst = stack(floods1, floods2, floods3, floods4, floods5, kmeans_floodmask)
plot(floodst)

Question 6:

plot_rsum = sum(floodst)

raster_stats = cellStats(floodst, stat=sum)*res(floodst)^2/1e6 #showing area covered in km&2

#Step 1

kable(raster_stats, col.names = c("Area Covered"), caption = "Area Covered by each classification") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Area Covered by each classification
Area Covered
layer.1 5.9994
layer.2 6.4908
layer.3 10.7451
layer.4 7.6221
layer.5 13.6809
layer.6 7.9371

Step 2

plot(plot_rsum, col = blues9)

# Step 3

plot_copy = (plot_rsum[plot_rsum == 0] = NA)
plot(plot_rsum, col = blues9)

Why is our raster STACK, comprised of 6 RASTERS, have values ranging from 6 to 0?

Extra Credit:

Used mapview from previous question to identify area of interest

aoi = st_point(c(-91.78946, 42.06303)) %>% 
  st_sfc(crs = 4326) %>% 
  st_sf() %>% 
  st_transform(st_crs(floodst))

aoi2 = st_point(c(-91.78482, 42.06245)) %>% 
  st_sfc(crs = 4326) %>% 
  st_sf() %>% 
  st_transform(st_crs(floodst))
# 
# extract(floodst, aoi)
# extract(floodst, aoi2)